Wolbachia (oligotyping)

Load packages, paths, functions

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("RColorBrewer", "ggpubr", "cowplot", "Biostrings", "openxlsx", "kableExtra")
invisible(lapply(packages, require, character.only = TRUE))

# create output folders if needed
if (!dir.exists("../../../../output/2_Oligotyping/2D")) {dir.create("../../../../output/2_Oligotyping/2D")}

Preparation

Tables preparation

Seqtab

# move to oligotyping directory
path_wolbachia <- "../../../../output/2_Oligotyping/2A/Wolbachia/2A_oligotyping_Wolbachia_sequences-c2-s1-a0.0-A0-M10"

# load the matrix count table
matrix_count <- read.table(paste0(path_wolbachia, "/MATRIX-COUNT.txt"), header = TRUE) %>% t()

# arrange it
colnames(matrix_count) <- matrix_count[1,]
matrix_count <- matrix_count[-1,]
matrix_count <- matrix_count %>% as.data.frame()

# print it
matrix_count %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP14 NP2 NP20 NP27 NP29 NP30 NP34 NP35 NP37 NP38 NP39 NP41 NP42 NP43 NP44 NP5 NP8 S100 S101 S102 S104 S105 S106 S107 S108 S109 S110 S113 S121 S122 S123 S124 S126 S127 S128 S146 S147 S148 S150 S151 S152 S153 S154 S160 S162 S163 S164 S165 S166 S167 S169 S170 S18 S19 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 S41 S42 S43 S44 S45 S46 S47 S48 S49 S50 S51 S52 S53 S55 S56 S57 S58 S59 S60 S61 S63 S64 S65 S79 S80 S83 S84 S85 S86 S87 S88 S89
AT 75 88 16 37 149 164 211 23 9 24 29 2 23 7510 618873 60 1132 9 64 29 637 101 5924 479 3 80 1571 2 703792 312332 52168 4204 3446 52042 55188 32800 51828 55473 58672 57141 11538 16105 7915 16231 14357 12030 19909 13241 20256 20238 11212 12053 18527 12233 13371 7439 57333 20325 9679 18121 14055 10723 10940 9242 7090 4234 30679 39952 22351 9576 35749 38103 47217 5222 5321 4331 18 4200 4782 2624 12821 5495 5240 4692 764 45671 11815 19172 7 4091 9243 7837 17710 498 1850 56207 33144 28127 60995 21063 4 50226 42384 48894 30175 32592 44266 49111 53121 47253 11 59372 52605 41938 56299 41407 61573 65420 52538 4
GT 10 15 4 10 25 41 42 5 2 1 7 0 3 195 15502 0 37 1 2 0 6 1 147 16 0 2 57 0 18486 14077 100 6 2 140 113 88 105 75 129 89 32 3040 1385 2732 2582 2206 3493 2227 3515 3517 2058 2253 3254 2047 2700 1570 10526 3468 1829 3061 2613 2266 2164 1511 1157 6 56 90 35 25 66 87 106 6 10 4 1 2 7 2 2 3 4 8 3 67 21 27 0 9 7 5 36 2 1 39 28 21 129 36 0 88 110 112 71 71 72 117 142 176 0 121 58 129 136 98 130 172 150 0
AG 2 1 1 1 4 9 6 0 5 2 2 0 2 186 11074 4 36 0 0 1 10 2 88 11 1 3 53 0 11433 6724 85 5 3 51 68 64 97 71 95 54 17 502 173 408 435 372 472 285 438 491 314 297 402 245 369 255 1623 592 314 431 454 450 367 352 254 6 35 31 31 12 66 45 55 8 8 4 0 8 10 3 12 7 8 2 2 38 17 22 0 1 5 5 7 0 3 60 39 23 67 29 0 33 25 53 29 46 37 42 32 70 0 49 42 56 63 41 41 69 67 0
TT 4 3 1 0 6 3 6 0 3 0 1 0 1 1 318 0 2 0 0 0 0 0 4 0 0 0 0 0 312 176 6 0 0 1 7 6 8 10 11 6 2 582 288 626 644 527 642 480 678 794 492 409 653 478 522 284 2539 660 418 620 476 420 448 369 292 0 3 7 2 2 7 4 7 0 2 0 0 0 2 0 2 1 1 0 0 4 2 2 0 1 4 1 2 0 0 3 1 4 8 2 0 6 2 7 2 5 6 4 13 9 0 10 3 7 12 8 3 8 4 0
CC 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1080 0 0 0 0 0 3 0 5 1 0 0 7 0 1082 608 73 12 3 116 108 68 87 84 93 72 21 18 9 47 42 26 48 18 30 31 18 16 22 24 18 18 130 31 15 38 31 22 19 12 16 5 63 84 31 17 51 99 109 8 9 4 0 3 7 3 19 4 3 3 3 62 21 32 0 4 20 12 41 0 2 69 35 16 125 52 0 75 81 74 58 57 78 236 97 83 0 138 74 87 113 76 152 177 170 0
GG 1 0 0 0 0 1 2 0 0 1 0 0 0 9 681 0 1 0 0 0 3 0 6 0 0 0 3 0 804 646 0 0 0 0 1 0 0 0 0 0 0 99 31 74 72 68 114 47 71 80 55 49 78 47 55 49 310 91 76 76 85 82 97 65 43 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0
AC 0 0 0 0 0 0 0 0 0 0 1 0 0 0 411 0 0 0 0 0 0 0 2 0 0 0 3 0 248 216 54 0 2 49 92 27 27 22 22 14 5 15 2 11 12 2 18 7 8 18 6 4 7 5 3 5 37 11 1 8 8 8 8 1 0 0 19 28 6 10 24 35 62 6 5 1 0 1 2 0 0 0 1 5 0 23 1 18 0 0 0 2 21 0 0 13 2 2 73 13 0 19 7 14 21 27 25 35 39 37 0 64 5 55 53 60 84 111 75 0

Taxonomy

# load the fasta table
fasta <- readDNAStringSet(paste0(path_wolbachia, "/OLIGO-REPRESENTATIVES.fasta"))

# arrange it
fasta <- fasta %>% as.data.frame()
colnames(fasta) <- "seq"
fasta$oligotype <- rownames(fasta)
fasta <- fasta %>% dplyr::select(-c(seq))

# print it
fasta %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
oligotype
AT AT
GT GT
AG AG
TT TT
CC CC
GG GG
AC AC

Change oligotype name by oligotype / MED nodes in the matrix count

# Reference file 

## load the reference table
ref_oligo_med2 <- read.table("../../../../output/2_Oligotyping/2B/2B_REF_info_wolbachia.tsv", sep="\t", header = TRUE)

## select only the 7 oligotypes of Wolbachia
ref_oligo_med2 <- ref_oligo_med2[!is.na(ref_oligo_med2$oligotype),]

## change order of columns
ref_oligo_med2 <- ref_oligo_med2 %>% select(c(seq, oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## create a column with reference name (will be used in plots)
ref_oligo_med2$ref <- paste0("oligotype_", ref_oligo_med2$OLIGO_oligotype_frequency_size, " / node_", ref_oligo_med2$MED_node_frequency_size)

## create a copy of fasta 
fasta2 <- fasta

# Matrix count

## create an oligotype column in the matrix count
matrix_count$oligotype <- rownames(matrix_count)

## change order of columns
matrix_count <- matrix_count %>% dplyr::select(c(oligotype, everything()))

## merge the matrix count and the reference dataframe
matrix_count2 <- matrix_count %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(c(oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size, ref, everything()))

## change rownames
rownames(matrix_count2) <- matrix_count2$ref

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(-c(oligotype, ref, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## print it
matrix_count2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP14 NP2 NP20 NP27 NP29 NP30 NP34 NP35 NP37 NP38 NP39 NP41 NP42 NP43 NP44 NP5 NP8 S100 S101 S102 S104 S105 S106 S107 S108 S109 S110 S113 S121 S122 S123 S124 S126 S127 S128 S146 S147 S148 S150 S151 S152 S153 S154 S160 S162 S163 S164 S165 S166 S167 S169 S170 S18 S19 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 S41 S42 S43 S44 S45 S46 S47 S48 S49 S50 S51 S52 S53 S55 S56 S57 S58 S59 S60 S61 S63 S64 S65 S79 S80 S83 S84 S85 S86 S87 S88 S89
oligotype_AC (80) | size:2504 / node_N0253 (80) | size:2481 0 0 0 0 0 0 0 0 0 0 1 0 0 0 411 0 0 0 0 0 0 0 2 0 0 0 3 0 248 216 54 0 2 49 92 27 27 22 22 14 5 15 2 11 12 2 18 7 8 18 6 4 7 5 3 5 37 11 1 8 8 8 8 1 0 0 19 28 6 10 24 35 62 6 5 1 0 1 2 0 0 0 1 5 0 23 1 18 0 0 0 2 21 0 0 13 2 2 73 13 0 19 7 14 21 27 25 35 39 37 0 64 5 55 53 60 84 111 75 0
oligotype_AG (109) | size:42030 / node_N0245 (109) | size:41133 2 1 1 1 4 9 6 0 5 2 2 0 2 186 11074 4 36 0 0 1 10 2 88 11 1 3 53 0 11433 6724 85 5 3 51 68 64 97 71 95 54 17 502 173 408 435 372 472 285 438 491 314 297 402 245 369 255 1623 592 314 431 454 450 367 352 254 6 35 31 31 12 66 45 55 8 8 4 0 8 10 3 12 7 8 2 2 38 17 22 0 1 5 5 7 0 3 60 39 23 67 29 0 33 25 53 29 46 37 42 32 70 0 49 42 56 63 41 41 69 67 0
oligotype_AT (120) | size:3890567 / node_N0711 (120) | size:3744518 75 88 16 37 149 164 211 23 9 24 29 2 23 7510 618873 60 1132 9 64 29 637 101 5924 479 3 80 1571 2 703792 312332 52168 4204 3446 52042 55188 32800 51828 55473 58672 57141 11538 16105 7915 16231 14357 12030 19909 13241 20256 20238 11212 12053 18527 12233 13371 7439 57333 20325 9679 18121 14055 10723 10940 9242 7090 4234 30679 39952 22351 9576 35749 38103 47217 5222 5321 4331 18 4200 4782 2624 12821 5495 5240 4692 764 45671 11815 19172 7 4091 9243 7837 17710 498 1850 56207 33144 28127 60995 21063 4 50226 42384 48894 30175 32592 44266 49111 53121 47253 11 59372 52605 41938 56299 41407 61573 65420 52538 4
oligotype_CC (92) | size:7065 / node_N0026 (92) | size:7069 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1080 0 0 0 0 0 3 0 5 1 0 0 7 0 1082 608 73 12 3 116 108 68 87 84 93 72 21 18 9 47 42 26 48 18 30 31 18 16 22 24 18 18 130 31 15 38 31 22 19 12 16 5 63 84 31 17 51 99 109 8 9 4 0 3 7 3 19 4 3 3 3 62 21 32 0 4 20 12 41 0 2 69 35 16 125 52 0 75 81 74 58 57 78 236 97 83 0 138 74 87 113 76 152 177 170 0
oligotype_GG (44) | size:4080 / node_N0250 (44) | size:4011 1 0 0 0 0 1 2 0 0 1 0 0 0 9 681 0 1 0 0 0 3 0 6 0 0 0 3 0 804 646 0 0 0 0 1 0 0 0 0 0 0 99 31 74 72 68 114 47 71 80 55 49 78 47 55 49 310 91 76 76 85 82 97 65 43 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0
oligotype_GT (111) | size:119651 / node_N0241 (111) | size:114576 10 15 4 10 25 41 42 5 2 1 7 0 3 195 15502 0 37 1 2 0 6 1 147 16 0 2 57 0 18486 14077 100 6 2 140 113 88 105 75 129 89 32 3040 1385 2732 2582 2206 3493 2227 3515 3517 2058 2253 3254 2047 2700 1570 10526 3468 1829 3061 2613 2266 2164 1511 1157 6 56 90 35 25 66 87 106 6 10 4 1 2 7 2 2 3 4 8 3 67 21 27 0 9 7 5 36 2 1 39 28 21 129 36 0 88 110 112 71 71 72 117 142 176 0 121 58 129 136 98 130 172 150 0
oligotype_TT (89) | size:15422 / node_N0244 (89) | size:15236 4 3 1 0 6 3 6 0 3 0 1 0 1 1 318 0 2 0 0 0 0 0 4 0 0 0 0 0 312 176 6 0 0 1 7 6 8 10 11 6 2 582 288 626 644 527 642 480 678 794 492 409 653 478 522 284 2539 660 418 620 476 420 448 369 292 0 3 7 2 2 7 4 7 0 2 0 0 0 2 0 2 1 1 0 0 4 2 2 0 1 4 1 2 0 0 3 1 4 8 2 0 6 2 7 2 5 6 4 13 9 0 10 3 7 12 8 3 8 4 0 
## edit the fasta dataframe
fasta2 <- fasta2 %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")
rownames(fasta2) <- fasta2$ref
fasta2 <- fasta2 %>% dplyr::select(-c(MED_node_frequency_size, OLIGO_oligotype_frequency_size, oligotype))

## print it
fasta2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
ref
oligotype_AC (80) | size:2504 / node_N0253 (80) | size:2481 oligotype_AC (80) | size:2504 / node_N0253 (80) | size:2481
oligotype_AG (109) | size:42030 / node_N0245 (109) | size:41133 oligotype_AG (109) | size:42030 / node_N0245 (109) | size:41133
oligotype_AT (120) | size:3890567 / node_N0711 (120) | size:3744518 oligotype_AT (120) | size:3890567 / node_N0711 (120) | size:3744518
oligotype_CC (92) | size:7065 / node_N0026 (92) | size:7069 oligotype_CC (92) | size:7065 / node_N0026 (92) | size:7069
oligotype_GG (44) | size:4080 / node_N0250 (44) | size:4011 oligotype_GG (44) | size:4080 / node_N0250 (44) | size:4011
oligotype_GT (111) | size:119651 / node_N0241 (111) | size:114576 oligotype_GT (111) | size:119651 / node_N0241 (111) | size:114576
oligotype_TT (89) | size:15422 / node_N0244 (89) | size:15236 oligotype_TT (89) | size:15422 / node_N0244 (89) | size:15236

Metadata

metadata <- read.csv("../../../../metadata/metadata.csv", sep=";")
rownames(metadata) <- metadata$Sample

Phyloseq object with oligotypes

# convert matrix_count into matrix and numeric
matrix_count <- matrix_count2 %>% as.matrix()
class(matrix_count) <- "numeric"

# phyloseq elements
OTU = otu_table(as.matrix(matrix_count), taxa_are_rows =TRUE)
TAX = tax_table(as.matrix(fasta2))
SAM = sample_data(metadata)

# phyloseq object
ps <- phyloseq(OTU, TAX, SAM)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 7 taxa and 120 samples ]
## sample_data() Sample Data:       [ 120 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 7 taxa by 1 taxonomic ranks ]
compute_read_counts(ps)
## [1] 4081319
# remove blanks
ps <- subset_samples(ps, Strain!="Blank")
ps <- check_ps(ps)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 7 taxa and 113 samples ]
## sample_data() Sample Data:       [ 113 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 7 taxa by 1 taxonomic ranks ]

Create new metadata with Percent

Load ps with all samples (for final plot)

ps.filter <- readRDS("../../../../output/1_MED/1D/1D_MED_phyloseq_decontam.rds")
ps.filter <- check_ps(ps.filter)

Edit new metadata with Percent_wolbachia

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0)))

# add read depth in sample table of phyloseq object
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)

# select Wolbachia
ps.wolbachia <- ps.filter %>% subset_taxa(Genus=="Wolbachia")

# add read depth of Wolbachia
sample_data(ps.filter)$Read_wolbachia <- sample_sums(ps.wolbachia)
sample_data(ps.filter) %>% colnames()
##  [1] "Sample"         "Strain"         "Field"          "Country"       
##  [5] "Organ"          "Species"        "Run"            "Control"       
##  [9] "Species_italic" "Strain_italic"  "Read_depth"     "is.neg"        
## [13] "Read_wolbachia"
sample_data(ps.wolbachia) %>% colnames()
##  [1] "Sample"         "Strain"         "Field"          "Country"       
##  [5] "Organ"          "Species"        "Run"            "Control"       
##  [9] "Species_italic" "Strain_italic"  "Read_depth"     "is.neg"
# add percent of Wolbachia
sample_data(ps.filter)$Percent_wolbachia <- sample_data(ps.filter)$Read_wolbachia / sample_data(ps.filter)$Read_depth

# round the percent of Wolbachia at 2 decimals
sample_data(ps.filter)$Percent_wolbachia <- sample_data(ps.filter)$Percent_wolbachia %>% round(2)

# extract metadata table
test <- data.frame(sample_data(ps.filter))

# merge this metadata table with the other
new.metadata <- data.frame(sample_data(ps)) %>% merge(test %>% dplyr::select(c(Sample, Read_depth, Read_wolbachia, Percent_wolbachia)), by="Sample")
new.metadata <- test[new.metadata$Sample %in% sample_names(ps),]
rownames(new.metadata) <- new.metadata$Sample

# print it
new.metadata %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
Sample Strain Field Country Organ Species Run Control Species_italic Strain_italic Read_depth is.neg Read_wolbachia Percent_wolbachia
CTC1 CTC1 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 29779 FALSE 0 0.00
CTC10 CTC10 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2609 FALSE 92 0.04
CTC11 CTC11 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 13874 FALSE 107 0.01
CTC12 CTC12 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1146 FALSE 22 0.02
CTC13 CTC13 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 18035 FALSE 48 0.00
CTC14 CTC14 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1708 FALSE 184 0.11
CTC15 CTC15 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 23180 FALSE 218 0.01
CTC2 CTC2 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 30692 FALSE 268 0.01
CTC3 CTC3 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 39920 FALSE 28 0.00
CTC4 CTC4 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2139 FALSE 19 0.01
CTC5 CTC5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 15789 FALSE 28 0.00
CTC6 CTC6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 19753 FALSE 40 0.00
CTC9 CTC9 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 4980 FALSE 29 0.01
NP14 NP14 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 7973 FALSE 7901 0.99
NP2 NP2 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 648335 FALSE 647941 1.00
NP20 NP20 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 136 FALSE 64 0.47
NP27 NP27 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 1234 FALSE 1208 0.98
NP29 NP29 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 203 FALSE 10 0.05
NP30 NP30 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 228 FALSE 66 0.29
NP34 NP34 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 95 FALSE 30 0.32
NP35 NP35 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 196532 FALSE 659 0.00
NP36 NP36 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 249 FALSE 0 0.00
NP37 NP37 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 419340 FALSE 104 0.00
NP38 NP38 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 282479 FALSE 6176 0.02
NP39 NP39 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 218684 FALSE 507 0.00
NP41 NP41 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 247152 FALSE 4 0.00
NP42 NP42 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 185157 FALSE 85 0.00
NP43 NP43 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 239335 FALSE 1694 0.01
NP44 NP44 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 156879 FALSE 2 0.00
NP5 NP5 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 736159 FALSE 736157 1.00
NP8 NP8 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 334799 FALSE 334782 1.00
S100 S100 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52486 FALSE 52486 1.00
S102 S102 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 3456 FALSE 3456 1.00
S104 S104 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52403 FALSE 52399 1.00
S105 S105 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 55577 FALSE 55577 1.00
S106 S106 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 33053 FALSE 33053 1.00
S107 S107 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52154 FALSE 52153 1.00
S108 S108 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 55735 FALSE 55735 1.00
S109 S109 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 59023 FALSE 59023 1.00
S110 S110 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 57377 FALSE 57377 1.00
S121 S121 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 20361 FALSE 20361 1.00
S122 S122 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 9803 FALSE 9803 1.00
S123 S123 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 20130 FALSE 20130 1.00
S124 S124 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 18146 FALSE 18145 1.00
S126 S126 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 15235 FALSE 15231 1.00
S127 S127 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 24696 FALSE 24696 1.00
S128 S128 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 16305 FALSE 16305 1.00
S146 S146 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 25012 FALSE 24996 1.00
S147 S147 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 25171 FALSE 25169 1.00
S148 S148 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 14164 FALSE 14155 1.00
S150 S150 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 15081 FALSE 15081 1.00
S151 S151 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 22944 FALSE 22944 1.00
S152 S152 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 15082 FALSE 15080 1.00
S153 S153 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 17040 FALSE 17038 1.00
S154 S154 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 9626 FALSE 9620 1.00
S160 S160 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 72508 FALSE 72501 1.00
S162 S162 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 25180 FALSE 25179 1.00
S163 S163 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 12333 FALSE 12332 1.00
S164 S164 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 22368 FALSE 22355 1.00
S165 S165 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 17731 FALSE 17723 1.00
S166 S166 Field - Camping Europe Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Camping~Europe 13979 FALSE 13971 1.00
S167 S167 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 14048 FALSE 14045 1.00
S169 S169 Field - Camping Europe Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Camping~Europe 11553 FALSE 11553 1.00
S170 S170 Field - Camping Europe Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Camping~Europe 8852 FALSE 8852 1.00
S18 S18 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 4290 FALSE 4251 0.99
S19 S19 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 44527 FALSE 30855 0.69
S20 S20 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 42864 FALSE 40192 0.94
S21 S21 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 33798 FALSE 22456 0.66
S22 S22 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 19044 FALSE 9642 0.51
S23 S23 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 38172 FALSE 35964 0.94
S24 S24 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 42355 FALSE 38373 0.91
S25 S25 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 47688 FALSE 47556 1.00
S26 S26 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 5394 FALSE 5250 0.97
S27 S27 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 24558 FALSE 5355 0.22
S28 S28 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 4503 FALSE 4344 0.96
S30 S30 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 25353 FALSE 4214 0.17
S31 S31 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 20417 FALSE 4810 0.24
S32 S32 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 12441 FALSE 2632 0.21
S33 S33 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 33867 FALSE 12856 0.38
S34 S34 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 9367 FALSE 5510 0.59
S35 S35 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 11663 FALSE 5257 0.45
S36 S36 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 33020 FALSE 4710 0.14
S37 S37 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 18340 FALSE 772 0.04
S38 S38 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 54790 FALSE 45865 0.84
S39 S39 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 36273 FALSE 11878 0.33
S40 S40 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 44448 FALSE 19274 0.43
S42 S42 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 4107 FALSE 4106 1.00
S43 S43 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 9279 FALSE 9279 1.00
S44 S44 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 8026 FALSE 7862 0.98
S45 S45 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 18150 FALSE 17817 0.98
S47 S47 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 1951 FALSE 1856 0.95
S48 S48 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 56738 FALSE 56391 0.99
S49 S49 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 33498 FALSE 33249 0.99
S50 S50 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 28481 FALSE 28193 0.99
S51 S51 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 61788 FALSE 61398 0.99
S52 S52 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 21553 FALSE 21195 0.98
S55 S55 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 50447 FALSE 50447 1.00
S56 S56 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 42609 FALSE 42609 1.00
S57 S57 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 49157 FALSE 49154 1.00
S58 S58 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 30357 FALSE 30357 1.00
S59 S59 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 32798 FALSE 32798 1.00
S60 S60 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 44485 FALSE 44485 1.00
S61 S61 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 49545 FALSE 49545 1.00
S63 S63 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 53444 FALSE 53444 1.00
S64 S64 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 47628 FALSE 47628 1.00
S79 S79 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 59755 FALSE 59755 1.00
S80 S80 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52788 FALSE 52788 1.00
S83 S83 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 42272 FALSE 42272 1.00
S84 S84 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 56676 FALSE 56676 1.00
S85 S85 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 41690 FALSE 41690 1.00
S86 S86 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 61984 FALSE 61984 1.00
S87 S87 Field - Bosc Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 65958 FALSE 65958 1.00
S88 S88 Field - Bosc Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 53102 FALSE 53004 1.00 
# replace metadata in the created phyloseq object
sample_data(ps) <- sample_data(new.metadata)

Taxonomic structure

Count

col <- brewer.pal(7, "Pastel2")

# reshape data for plot
test3 <- test %>% select(c(Sample, Species, Strain, Organ, Read_depth, Read_wolbachia)) %>% reshape2::melt(id.vars=c("Sample", "Species", "Strain", "Organ"), vars=c("Read_depth", "Read_wolbachia"))

count_whole <- test3[test3$Organ=="Whole",]
count_ovary <- test3[test3$Organ=="Ovary",]

make.italic <- function(x) as.expression(lapply(x, function(y) bquote(italic(.(y)))))

levels(count_whole$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_ovary$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_whole$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

levels(count_ovary$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))


# plot
p_count1 <- ggplot(count_whole, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=12, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=1.1), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
p_count2 <- ggplot(count_ovary, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
 facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
    ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=0.8), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
# afficher plot
p_count1
## Warning: position_dodge requires non-overlapping x intervals
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

p_count2

# panels
p_group <- plot_grid(p_count1+theme(legend.position="none"), 
          p_count2+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("B1", "B2"), c(0, 0), c(1, 0.5), size = 20)
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
legend_plot <- get_legend(p_count1 + theme(legend.position="bottom"))
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_counts

Whole (the most abundant nodes)

#fig.height = 12, fig.width=18
guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0),
                                            nrow=4, byrow=TRUE))

# select whole
ps.filter.whole <- subset_samples(ps, Organ=="Whole")
ps.filter.whole <- prune_taxa(taxa_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole <- prune_samples(sample_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 7 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 14 sample variables ]
## tax_table()   Taxonomy Table:    [ 7 taxa by 1 taxonomic ranks ]
# data pour plot
#data_for_plot2 <- taxo_data_fast(ps.filter.whole, method = "abundance")
data_for_plot2 <- taxo_data(ps.filter.whole, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot2$Name), "\",\n") %>% cat()
## "oligotype_AC (80) | size:2504 / node_N0253 (80) | size:2481.",
##  "oligotype_AG (109) | size:42030 / node_N0245 (109) | size:41133.",
##  "oligotype_AT (120) | size:3890567 / node_N0711 (120) | size:3744518.",
##  "oligotype_CC (92) | size:7065 / node_N0026 (92) | size:7069.",
##  "oligotype_GG (44) | size:4080 / node_N0250 (44) | size:4011.",
##  "oligotype_GT (111) | size:119651 / node_N0241 (111) | size:114576.",
##  "oligotype_TT (89) | size:15422 / node_N0244 (89) | size:15236.",
data_for_plot2$Name <- data_for_plot2$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot2$Name <- as.factor(data_for_plot2$Name)

new_names <- c( "oligotype_AT (120) | size:3890567 / N0711 (120) | size:3744518.",
                "oligotype_GT (111) | size:119651 / N0241 (111) | size:114576.",
                "oligotype_AG (109) | size:42030 / N0245 (109) | size:41133.",
                "oligotype_TT (89) | size:15422 / N0244 (89) | size:15236.",
                "oligotype_CC (92) | size:7065 / N0026 (92) | size:7069.",
                "oligotype_GG (44) | size:4080 / N0250 (44) | size:4011.",
                "oligotype_AC (80) | size:2504 / N0253 (80) | size:2481."
)

data_for_plot2$Name <- factor(data_for_plot2$Name, levels = new_names)

col <- c("oligotype_AT (120) | size:3890567 / N0711 (120) | size:3744518."="#FEB24C", 
         "oligotype_GT (111) | size:119651 / N0241 (111) | size:114576."="#FAD769", 
         "oligotype_AG (109) | size:42030 / N0245 (109) | size:41133."="#666666", 
         "oligotype_CC (92) | size:7065 / N0026 (92) | size:7069."="#BEAED4", 
         "oligotype_TT (89) | size:15422 / N0244 (89) | size:15236."="#BF5B17",
         "oligotype_AC (80) | size:2504 / N0253 (80) | size:2481."="#F4CAE4", 
         "oligotype_GG (44) | size:4080 / N0250 (44) | size:4011."="#FDCDAC")

levels(data_for_plot2$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(data_for_plot2$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

#data_for_plot2 <- data_for_plot2 %>% na.omit()

p2 <- ggplot(data_for_plot2, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p2

Ovary (the most abundant nodes)

# select ovary
ps.filter.ovary <- subset_samples(ps, Organ=="Ovary")
ps.filter.ovary <- prune_taxa(taxa_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary <- prune_samples(sample_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 7 taxa and 46 samples ]
## sample_data() Sample Data:       [ 46 samples by 14 sample variables ]
## tax_table()   Taxonomy Table:    [ 7 taxa by 1 taxonomic ranks ]
# data pour plot
#data_for_plot3 <- taxo_data_fast(ps.filter.ovary, method = "abundance")
data_for_plot3 <- taxo_data(ps.filter.ovary, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot3$Name), "\",\n") %>% cat()
## "oligotype_AC (80) | size:2504 / node_N0253 (80) | size:2481.",
##  "oligotype_AG (109) | size:42030 / node_N0245 (109) | size:41133.",
##  "oligotype_AT (120) | size:3890567 / node_N0711 (120) | size:3744518.",
##  "oligotype_CC (92) | size:7065 / node_N0026 (92) | size:7069.",
##  "oligotype_GG (44) | size:4080 / node_N0250 (44) | size:4011.",
##  "oligotype_GT (111) | size:119651 / node_N0241 (111) | size:114576.",
##  "oligotype_TT (89) | size:15422 / node_N0244 (89) | size:15236.",
data_for_plot3$Name <- data_for_plot3$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot3$Name <- as.factor(data_for_plot3$Name)

new_names <- c( "oligotype_AT (120) | size:3890567 / N0711 (120) | size:3744518.",
                "oligotype_GT (111) | size:119651 / N0241 (111) | size:114576.",
                "oligotype_AG (109) | size:42030 / N0245 (109) | size:41133.",
                "oligotype_TT (89) | size:15422 / N0244 (89) | size:15236.",
                "oligotype_CC (92) | size:7065 / N0026 (92) | size:7069.",
                "oligotype_GG (44) | size:4080 / N0250 (44) | size:4011.",
                "oligotype_AC (80) | size:2504 / N0253 (80) | size:2481."
)

data_for_plot3$Name <- factor(data_for_plot3$Name, levels = new_names)

levels(data_for_plot3$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

#levels(data_for_plot3$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))
levels(data_for_plot3$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)")

#data_for_plot3 <- data_for_plot3 %>% na.omit()

p3 <- ggplot(data_for_plot3, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller = label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p3

Panels taxonomy of whole / ovary

legend_plot <- get_legend(p2 + theme(legend.position="bottom"))

# panels
p_group <- plot_grid(p2+theme(legend.position="none"), 
          p3+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("A1", "A2"), c(0, 0), c(1, 0.5), size = 20)

p_taxo <- plot_grid(p_group, legend_plot, nrow=2, rel_heights = c(1, .1))
p_taxo

Save taxonomic plot

tiff("../../../../output/2_Oligotyping/2D/2Da_OLIGO_counts_wolbachia.tiff", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/2_Oligotyping/2D/2Da_OLIGO_taxonomic_wolbachia_whole.tiff", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/2_Oligotyping/2D/2Da_OLIGO_taxonomic_wolbachia_ovary.tiff", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/2_Oligotyping/2D/2Da_OLIGO_taxonomic_wolbachia.tiff", units="in", width=18, height=16, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Da_OLIGO_counts_wolbachia_big.png", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Da_OLIGO_counts_wolbachia_small.png", units="in", width=18, height=14, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Da_OLIGO_taxonomic_wolbachia_whole.png", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Da_OLIGO_taxonomic_wolbachia_ovary.png", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Da_OLIGO_taxonomic_wolbachia_big.png", units="in", width=18, height=18, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Da_OLIGO_taxonomic_wolbachia_small.png", units="in", width=18, height=14, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2

Make main plot

img <- magick::image_read(paste0(path_wolbachia, "/HTML-OUTPUT/entropy.png"))
p_entropy <- magick::image_ggplot(img, interpolate = TRUE)
p_entropy+ theme(plot.margin = unit(c(-7,-2.5,-7,-0.5), "cm"))

p_entropy+ theme(plot.margin=unit(c(-7,-2,-12,-5), "mm"))

aligned <- plot_grid(p_taxo, 
                     p_counts, 
                     align="hv")

aligned

p_entropy2 <- plot_grid(p_entropy, nrow=1)+
  draw_plot_label(c("C"), c(0), c(1), size=20, hjust=-0.5)

p_entropy2

t_plot <- plot_grid(aligned, 
                    p_entropy2,
                    nrow=2, 
                    ncol=1, 
                    scale=1,
                    rel_heights=c(2,1))

t_plot

tiff("../../../../output/2_Oligotyping/2D/2Da_OLIGO_main_wolbachia.tiff", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Da_OLIGO_main_wolbachia.png", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2

Session info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.3.1    openxlsx_4.2.3      Biostrings_2.54.0  
##  [4] XVector_0.26.0      IRanges_2.20.2      S4Vectors_0.24.4   
##  [7] BiocGenerics_0.32.0 cowplot_1.1.0       ggpubr_0.4.0       
## [10] RColorBrewer_1.1-2  forcats_0.5.0       stringr_1.4.0      
## [13] dplyr_1.0.2         purrr_0.3.4         readr_1.4.0        
## [16] tidyr_1.1.2         tibble_3.0.4        tidyverse_1.3.0    
## [19] ggplot2_3.3.2       phyloseq_1.30.0    
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-0  ggsignif_0.6.0    ellipsis_0.3.1    rio_0.5.16       
##  [5] fs_1.5.0          rstudioapi_0.11   farver_2.0.3      fansi_0.4.1      
##  [9] lubridate_1.7.9   xml2_1.3.2        codetools_0.2-16  splines_3.6.3    
## [13] knitr_1.30        ade4_1.7-15       jsonlite_1.7.1    broom_0.7.2      
## [17] cluster_2.1.0     dbplyr_1.4.4      compiler_3.6.3    httr_1.4.2       
## [21] backports_1.1.10  assertthat_0.2.1  Matrix_1.2-18     cli_2.1.0        
## [25] htmltools_0.5.1.1 tools_3.6.3       igraph_1.2.6      gtable_0.3.0     
## [29] glue_1.4.2        reshape2_1.4.4    Rcpp_1.0.5        carData_3.0-4    
## [33] Biobase_2.46.0    cellranger_1.1.0  jquerylib_0.1.3   vctrs_0.3.4      
## [37] multtest_2.42.0   ape_5.4-1         nlme_3.1-149      iterators_1.0.13 
## [41] xfun_0.22         ps_1.4.0          rvest_0.3.6       lifecycle_0.2.0  
## [45] rstatix_0.6.0     zlibbioc_1.32.0   MASS_7.3-53       scales_1.1.1     
## [49] hms_0.5.3         biomformat_1.14.0 rhdf5_2.30.1      yaml_2.2.1       
## [53] curl_4.3          sass_0.3.1        stringi_1.5.3     highr_0.8        
## [57] foreach_1.5.1     permute_0.9-5     zip_2.1.1         rlang_0.4.10     
## [61] pkgconfig_2.0.3   evaluate_0.14     lattice_0.20-41   Rhdf5lib_1.8.0   
## [65] labeling_0.4.2    tidyselect_1.1.0  plyr_1.8.6        magrittr_1.5     
## [69] bookdown_0.22     R6_2.4.1          magick_2.5.2      generics_0.0.2   
## [73] DBI_1.1.0         pillar_1.4.6      haven_2.3.1       foreign_0.8-75   
## [77] withr_2.3.0       mgcv_1.8-33       survival_3.2-7    abind_1.4-5      
## [81] modelr_0.1.8      crayon_1.3.4      car_3.0-10        rmarkdown_2.7    
## [85] grid_3.6.3        readxl_1.3.1      data.table_1.13.2 blob_1.2.1       
## [89] vegan_2.5-6       rmdformats_1.0.2  reprex_0.3.0      digest_0.6.26    
## [93] webshot_0.5.2     munsell_0.5.0     viridisLite_0.3.0 bslib_0.2.4